Causal vs Acausal - Slide from Openmodelica presentation
CirculatorySystemModels.jl
CirculatorySystemModels.jl (Schenkel et al.) is a composable modelling framework for cardiovascular circulation models using ModelingToolkit.jl.
It allows the straighforward generation of CVS models from a component library and automates the creation of ODE or ODAE systems for arbitrary flow networks.
%%{init: {'theme': 'base', 'themeVariables': { 'clusterBkg': '#ddd'}}}%%
flowchart
A[Network] --> C(Connections)
B[Components] --> D(Component Equations)
C --> E[System Equations]
D --> E
E --> |Structural Simplify| G[ODEProblem]
G --> H[ODESolver]
Another 2x speed-up compared to the hand-derived ODEs! MTK optimises ODEs for efficiency and stability. In this case ODEs are not for \(p_{1}\) and \(p_{2}\), but for \(Q_{L_p}\) and \(\Delta p_C\).
`"""`ShiSystemicLoop(; name, SAS_C, SAS_R, SAS_L, SAT_C, SAT_R, SAT_L, SAR_R, SCP_R, SVN_C, SVN_R)`Implements systemic loop as written by Shi in [Shi].Parameters are in the cm, g, s system.Pressure in mmHg.Volume in ml.Flow in cm^3/s (ml/s).Named parameters:`SAS_C`: Aortic sinus compliance in ml/mmHg`SAS_R`: Aortic sinus resistance in mmHg*s/ml`SAS_L`: Aortic sinus inductance in mmHg*s^2/ml`SAT_C`: Artery compliance in ml/mmHg`SAT_R`: Artery resistance in mmHg*s/ml`SAT_L`: Artery inductance in mmHg*s^2/ml`SAR_R`: Arteriole resistance in mmHg*s/ml`SCP_R`: Capillary resistance in mmHg*s/ml`SVN_C`: Vein compliance in ml/mmHg`SVN_R`: Vein resistance in mmHg*s/ml"""function ShiSystemicLoop(; name, SAS_C, SAS_R, SAS_L, SAT_C, SAT_R, SAT_L, SAR_R, SCP_R, SVN_C, SVN_R) @named in = Pin() @named out = Pin() sts = @variables Δp(t) = 0.0 q(t) = 0.0 ps = [] # No parameters in this function # These are the components the subsystem is made of: ## Systemic Aortic Sinus ## @named SAS = CRL(C=SAS_C, R=SAS_R, L=SAS_L) ## Systemic Artery ## @named SAT = CRL(C=SAT_C, R=SAT_R, L=SAT_L) ## Systemic Arteriole ## @named SAR = Resistor(R=SAR_R) ## Systemic Capillary ## @named SCP = Resistor(R=SCP_R) ## Systemic Vein ## @named SVN = CR(C=SVN_C, R=SVN_R) # The equations for the subsystem are created by # 'connect'-ing the components eqs = [ Δp ~ out.p - in.p q ~ in.q connect(in, SAS.in) connect(SAS.out, SAT.in) connect(SAT.out, SAR.in) connect(SAR.out, SCP.in) connect(SCP.out, SVN.in) connect(SVN.out, out) ] # and finaly compose the system compose(ODESystem(eqs, t, sts, ps; name=name), in, out, SAS, SAT, SAR, SCP, SVN)end
`ShiHeart(; name, τ, LV_V₀, LV_p0, LV_Emin, LV_Emax, LV_τes, LV_τed, LV_Eshift, RV_V₀, RV_p0, RV_Emin, RV_Emax, RV_τes, RV_τed, RV_Eshift, LA_V₀, LA_p0, LA_Emin, LA_Emax, LA_τes, LA_τed, LA_Eshift, RA_V₀, RA_p0, RA_Emin, RA_Emax, RA_τes, RA_τed, RA_Eshift, AV_CQ, AV_Kp, AV_Kf, AV_Kb, AV_Kv, AV_θmax, AV_θmin, PV_CQ, PV_Kp, PV_Kf, PV_Kb, PV_Kv, PV_θmax, PV_θmin, MV_CQ, MV_Kp, MV_Kf, MV_Kb, MV_Kv, MV_θmax, MV_θmin, TV_CQ, TV_Kp, TV_Kf, TV_Kb, TV_Kv, TV_θmax, TV_θmin)`Models a whole heart, made up of 2 ventricles (Left & Right Ventricle) and 2 atria (Left & Right atrium)created from the ShiChamber element. Includes the 4 corresponding valves (Aortic, Mitral, Pulmonary and Tricuspid valve) created using the ShiValve element.Parameters are in the cm, g, s system.Pressure in mmHg.Volume in ml.Flow in cm^3/s (ml/s).Maximum and Minimum angles given in rad, to convert from degrees multiply angle by pi/180.Named parameters:`τ` Length of the cardiac cycle in s`LV_V₀` Unstressed left ventricular volume in ml`LV_p0` Unstressed left ventricular pressure in mmHg`LV_Emin` Minimum left ventricular elastance (diastole) in mmHg/ml`LV_Emax` Maximum left ventricular elastance (systole) in mmHg/ml`LV_τes` Left ventricular end systolic time in s`LV_τed` Left ventricular end distolic time in s`LV_Eshift` Shift time of contraction - 0 for left ventricle`RV_V₀` Unstressed right ventricular volume in ml`RV_p0` Unstressed right ventricular pressure in mmHg`RV_Emin` Minimum right ventricular elastance (diastole) in mmHg/ml`RV_Emax` Maximum right ventricular elastance (systole) in mmHg/ml`RV_τes` Right ventricular end systolic time in s`RV_τed` Right ventricular end distolic time in s`RV_Eshift` Shift time of contraction - 0 for right ventricle`LA_V₀` Unstressed left atrial volume in ml`LA_p0` Unstressed left atrial pressure in mmHg`LA_Emin` Minimum left atrial elastance (diastole) in mmHg/ml`LA_Emax` Maximum left atrial elastance (systole) in mmHg/ml`LA_τes` Left atrial end systolic time in s`LA_τed` Left atrial end distolic time in s`LA_Eshift` Shift time of contraction in s`RA_V₀` Unstressed right atrial volume in ml`RA_p0` Unstressed right atrial pressure in mmHg`RA_Emin` Minimum right atrial elastance (diastole) in mmHg/ml`RA_Emax` Maximum right atrial elastance (systole) in mmHg/ml`RA_τes` Right atrial end systolic time in s`RA_τed` Right atrial end distolic time in s`RA_Eshift` Shift time of contraction in s`AV_CQ` Aortic valve flow coefficent in ml/(s*mmHg^0.5)`AV_Kp` Pressure effect on the aortic valve in rad/(s^2*mmHg)`AV_Kf` Frictional effect on the aortic valve in 1/s`AV_Kb` Fluid velocity effect on the aortic valve in rad/(s*m)`AV_Kv` Vortex effect on the aortic valve in rad/(s*m)`AV_θmax` Aortic valve maximum opening angle in rad`AV_θmin` Aortic valve minimum opening angle in rad`MV_CQ` Mitral valve flow coefficent in ml/(s*mmHg^0.5)`MV_Kp` Pressure effect on the mitral valve in rad/(s^2*mmHg)`MV_Kf` Frictional effect on the mitral valve in 1/s`MV_Kb` Fluid velocity effect on the mitral valve in rad/(s*m)`MV_Kv` Vortex effect on the mitral valve in rad/(s*m)`MV_θmax` Mitral valve maximum opening angle in rad`MV_θmin` Mitral valve minimum opening angle in rad`PV_CQ` Pulmonary valve flow coefficent in ml/(s*mmHg^0.5)`PV_Kp` Pressure effect on the pulmonary valve in rad/(s^2*mmHg)`PV_Kf` Frictional effect on the pulmonary valve in 1/s`PV_Kb` Fluid velocity effect on the pulmonary valve in rad/(s*m)`PV_Kv` Vortex effect on the pulmonary valve in rad/(s*m)`PV_θmax` Pulmonary valve maximum opening angle in rad`PV_θmin` Pulmonary valve minimum opening angle in rad`TV_CQ` Tricuspid valve flow coefficent in ml/(s*mmHg^0.5)`TV_Kp` Pressure effect on the tricuspid valve in rad/(s^2*mmHg)`TV_Kf` Frictional effect on the tricuspid valve in 1/s`TV_Kb` Fluid velocity effect on the tricuspid valve in rad/(s*m)`TV_Kv` Vortex effect on the pulmonary valve in rad/(s*m)`TV_θmax` Tricuspid valve maximum opening angle in rad`TV_θmin` Tricuspid valve minimum opening angle in rad"""function ShiHeart(; name, τ, LV_V₀, LV_p0, LV_Emin, LV_Emax, LV_τes, LV_τed, LV_Eshift, RV_V₀, RV_p0, RV_Emin, RV_Emax, RV_τes, RV_τed, RV_Eshift, LA_V₀, LA_p0, LA_Emin, LA_Emax, LA_τes, LA_τed, LA_Eshift, RA_V₀, RA_p0, RA_Emin, RA_Emax, RA_τes, RA_τed, RA_Eshift, AV_CQ, AV_Kp, AV_Kf, AV_Kb, AV_Kv, AV_θmax, AV_θmin, PV_CQ, PV_Kp, PV_Kf, PV_Kb, PV_Kv, PV_θmax, PV_θmin, MV_CQ, MV_Kp, MV_Kf, MV_Kb, MV_Kv, MV_θmax, MV_θmin, TV_CQ, TV_Kp, TV_Kf, TV_Kb, TV_Kv, TV_θmax, TV_θmin) @named LHin = Pin() @named LHout = Pin() @named RHin = Pin() @named RHout = Pin() sts = @variables Δp(t) = 0.0 q(t) = 0.0 # sts = [] # ps = @parameters Rc=Rc Rp=Rp C=C # No parameters in this function # Parameters are inherited from subcomponents ps = [] # These are the components the subsystem is made of: # Ventricles and atria @named LV = ShiChamber(V₀=LV_V₀, p₀=LV_p0, Eₘᵢₙ=LV_Emin, Eₘₐₓ=LV_Emax, τ=τ, τₑₛ=LV_τes, τₑₚ=LV_τed, Eshift=LV_Eshift) @named LA = ShiChamber(V₀=LA_V₀, p₀=LA_p0, Eₘᵢₙ=LA_Emin, Eₘₐₓ=LA_Emax, τ=τ, τₑₛ=LA_τes, τₑₚ=LA_τed, Eshift=LA_Eshift) @named RV = ShiChamber(V₀=RV_V₀, p₀=RV_p0, Eₘᵢₙ=RV_Emin, Eₘₐₓ=RV_Emax, τ=τ, τₑₛ=RV_τes, τₑₚ=RV_τed, Eshift=RV_Eshift) @named RA = ShiChamber(V₀=RA_V₀, p₀=RA_p0, Eₘᵢₙ=RA_Emin, Eₘₐₓ=RA_Emax, τ=τ, τₑₛ=RA_τes, τₑₚ=RA_τed, Eshift=RA_Eshift) # Valves @named AV = ShiValve(CQ=AV_CQ, Kp=AV_Kp, Kf=AV_Kf, Kb=AV_Kb, Kv=AV_Kv, θmax=AV_θmax, θmin=AV_θmin) @named MV = ShiValve(CQ=MV_CQ, Kp=MV_Kp, Kf=MV_Kf, Kb=MV_Kb, Kv=MV_Kv, θmax=MV_θmax, θmin=MV_θmin) @named TV = ShiValve(CQ=TV_CQ, Kp=TV_Kp, Kf=TV_Kf, Kb=TV_Kb, Kv=TV_Kv, θmax=TV_θmax, θmin=TV_θmin) @named PV = ShiValve(CQ=PV_CQ, Kp=PV_Kp, Kf=PV_Kf, Kb=PV_Kb, Kv=PV_Kv, θmax=PV_θmax, θmin=PV_θmin) # The equations for the subsystem are created by # 'connect'-ing the components eqs = [ Δp ~ out.p - in.p q ~ in.q connect(LHin, LA.in) connect(LA.out, MV.in) connect(MV.out, LV.in) connect(LV.out, AV.in) connect(AV.out, LHout) connect(RHin, RA.in) connect(RA.out, TV.in) connect(TV.out, RV.in) connect(RV.out, PV.in) connect(PV.out, RHout) ] # and finaly compose the system compose(ODESystem(eqs, t, sts, ps; name=name), in, out, LHin, LHout, RHin, RHout, LV, RV, LA, RA, AV, MV, TV, PV)end
Change valve model by changing 4 lines of code2. Introduce stenosis by changing stenosis parameter in 1 line of code. Quick simulation of variations (0.3s for Shi valve, 2.5s for new mechanistic, non-linear valve model).
Comparison to CellML
CellML is a modelling framework for biological modelling that has been developed since the early 1990s. While it is still being used, development seems to have stagnated somewhat. The main runtime engine for CellML is OpenCOR, which is implemented in C++, but limited to CellML v1.03.
I could not find a cardiovascular model that would run in OpenCOR, so for comparison of performance we compare to models that were tested in the Julia implementation of CellML CellMLToolkit: the simpler, electrophysiological models, O’Hara-Rudy, and ten Tusscher-Noble-Noble-Panfilov.
Model
Solver
OpenCOR
Julia
OR
CVODE
31ms
23ms
tTNNP
CVODE
7ms
8.3ms
So we have parity between those two implementations of the CellML modelling framework, one in C++, the other in Julia.
can read majority of CellML models and convert to ModelingToolkit (MTK) models5
equal performance as native C++ implementation in OpenCOR
integration into wider SciML framework: Optimisation, Identifiability, Global Sensitivity
very active developer community which are approachable
Key Disadvantages:
not yet established in the field
not unit aware (current development will change this!)
References
Korakianitis, Theodosios, and Yubing Shi. 2006. “A Concentrated Parameter Model for the Human Cardiovascular System Including Heart Valve Dynamics and Atrioventricular Interaction.”Medical Engineering & Physics 28 (7): 613–28.
Shi, Yubing, Patricia Lawford, and D. Rodney Hose. 2018. “Construction of Lumped-Parameter Cardiovascular Models Using the CellML Language.”Journal of Medical Engineering & Technology 42 (7): 525–31.
Footnotes
CellML does not support complex valve types.
Not including the definition of the valve mechanics itself.
Note: The C++ framework is not as flexible as the Julia framework, though, and many of the more advanced cellML models can only be run in Julia, since many of the tools in cellML seem to be stagnant and unsupported. OpenCell had it’s last release in 2010 and I could not get it to work on either Windows or Linux. OpenCOR is limited to cellML format 1.0 while most of the models in the repository are version 1.1.
at least in my opinion
In fact, the Julia version of CellML is the only one I could get to work on the cardiovascular system models! OpenCell is outdated and does not run on modern OSs, OpenCOR is limited to CellML v1.0.